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In the framework of the hydrodynamic model describing the response of electrons in a metal, we 
show that arrays of very narrow and shallow metallic slits have an optical response that is influenced 
by the spatial dispersion in metals arising from the repulsive interaction between electrons. As a 
simple Fabry-Perot model is not accurate enough to describe the structure’s behavior, we propose 
to consider the slits as generalized cavities with two modes, one being propagative and the other 
evanescent. This very general model allows to conclude that the impact of spatial dispersion on 
the propagative mode is the key factor explaining why the whole structure is sensitive to spatial 
dispersion. As the fabrication of such structures with relatively large gaps compared to previous 
experiments is within onr reach, this work paves the way for future much needed experiments on 
nonlocality. 


I. INTRODUCTION 

Drude’s modeP has proven unbelievably accurate 
throughout the twentieth century to describe the opti¬ 
cal response of metals, despite extensive studies in the 
seventies and eighties to find its limit ^2131 Recent ex¬ 
periments have however shown that the behavior of reso¬ 
nances in sub-nanometer-sized gaps can not be explained 
with Drude’s model alond^, making it necessary to take 
into account the repulsion between electrons inside the 
metal in the framework of a hydrodynamic modeP^^. In 
that case, the metallic response is spatially dispersive 
and since it cannot be reduced to a single permittivity 
depending on the frequency alone, it is often said to be 
non-local. Further experiments would however be wel¬ 
come because of th e sub-nanometer dimensions of the 
considered gap^^Hni, 

In the present work, we explain why the non-local re¬ 
sponse of metals can be expected to have an impact on 
deeply subwavelength metallic gratings, a structure that 
has been extensively stud ied in the past decade for its 
extraordinary transmissiorfSE^ and absorpt iorP^lli^. We 
underline that these effect will be clear for grooves that 
are as large as a few nanometers and that the fabrication 
of such still very narrow slits seems totally within our 

reachiSlli^. 

In the first part of this article, we will focus on the 
physical analysis of the absorption by the metallic slits ar¬ 
ray when spatial dispersion is neglected. It is now well ac¬ 
cepted that the extraordinary optical transmission of slit 
arrays i s due to the excitation of cavity resonances inside 
the slit^SE^, even if non-resonant mechanism allow for 
a high transmission for very thin structureJi^. The only 
guided mode propagating in the slits in p-polarization 
has actually no cut-off and can propagate whatever the 
slit width, with a wavevector kz close to fcp, the wavevec- 


tor in vacuum in most cases. This is why, as is usual for 
cavities, the thickness of the grating is roughly half the 
wavelength in vacuum for the fundamental resonanc^^, 
except for exotic case^^. The same mechanism explains 
the absorption by subwavelength grooves, the difference 
being that, since the cavity is now closed on one end (see 
Fig. [^, the depth of the grooves is only a quarter of the 
wavelength in vacuunfSl all cases, an extraordinar¬ 
ily strong tunneling effect expl ains th e way the energy is 
literally sucked up into the slit^SE^. 

As long as the m etal can be considered as almost per¬ 
fect (in the IRpE! or THJ^ or even microwav^^ range) 
or when the slit is larger than 50 nm in the optical range, 
the above physical picture is fully accurate. However, 
below that 50-nm threshold for the slit width in the op¬ 
tical range, the guided mode is more localized in the 
metal than in the slits because the skin depth S is typ¬ 
ically around 25 nm for noble metals. The mode then 
experiences what can be called a plasmonic drag: as 
the width of the gap decreases, it becomes slower and 
slower, its group velocity goes to zero and its wavevec¬ 
tor kz diverge^^. Its effective wavelength, defined as 
Ae// = 27:/kz thus becomes very small. The depth 
of the cavity being actually proportional to this effec¬ 
tive wavelengtlPSl, the thickness of the grating can be 
made extremely small and the slits will still constitute 
a resonatorfi^. This is the regime we are interested in, 
because it is one of the rare structures in which a plas¬ 
monic guided mode with a very high wavevector can be 
excited^. We show however that the physics of these 
resonators is slightly more complicated that previously 
thoughlP^ the slits are so shallow that the guided mode 
is not the only channel for light t o reac h the bottom of 
the structure. A one-mode modepdll^ is thus not suf¬ 
ficient to describe the cavity accurately. We give here 
a generalized Fabry-Perot formula to better predict the 
behavior of the resonances. 
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FIG. 1. (Color online) Representation of the grooves carved 
in metal, illuminated from above by a TM-polarized plane 
wave. We distinguish three different layers, layer II being the 
layer containing the slits. 


In a second part, we take into account spatial disper¬ 
sion in the framew ork o f the hydrodynamic model with 
realistic parameter^^^, and show that the response of 
the structure is influenced by nonlocality. Moreover, us¬ 
ing the generalized Fabry-Perot formula, we show that 
the influence of nonlocality on the wavevector of the 
guided mode explains almost totally why the structure 
is so sensitive to spatial dispersion in metals, completing 
the physical picture. 


II. GENERALIZED CAVITY MODEL 


The structure considered her j^^hs ue] jg presented in 
Figure It is an infinite array of deeply subwavelength 
grooves of width a ranging 2 to 5 nm and depth /i, from 
15 to 30 nm typically. The period d is ranging from 
30 to 50 nm in the present study. Since the structure 
is periodic with a period of the order of twice the skin 
depth, the slits must be considered as coupled, and the 
modes propagating in the slits can be found for normal 
incidence by solving the classical dispersion relation for 
metallo-dielectric structureJ^^ 


cosh(Ka) cosh{Kt{d — a)) 


Kt 

en 




sinh(Ka) sinh {Kt{d — a)) — 1 = 0 , ( 1 ) 


where — ek^ and = k"^ — kg, e being the 

permittivity of the metal and fcg = 27r/A. 

The different modes are indexed by n and characterized 
by a propagation constant kz along the z axis, and a 
magnetic field profile Hy{x). The magnetic field in layer 
II can thus be written as a sum of modes propagating 
upward or downward 


Hy{x,z) = Y,H'^{x) 6+*'="^] 

n 

( 2 ) 

As expected, there is only one propagating Bloch mode 
(presenting a propagation constant fci with a dominant 
real part). All the other modes are evanescent and thus 
attenuated in the 0 direction, with essentially imaginary 
propagation constants fc„. The propagative mode is re¬ 
flected back and forth in the slits, thus producing the 
resonance. The actual reflection coefficients (ri for the 
interface between the grating and air, r\ for th e bott om 
of the grooves) can be computed usi ng as 

is quite common for metallic gratingsl^^^U y \/0 under¬ 
line that the computation of the reflection coefficients re¬ 
quires the computation of many evanescent modes. One 
could expect, from the vast literature on the subject, that 
the reflection coefficient of the whole structure can simply 
be written using a Fabry-Perot formula: 


r = ro 


ri tpi tio 


( 3 ) 


where toi is the transmission coefficient from the incom¬ 
ing plane wave to the propagating mode in the slits and 
tio the transmission coefficient from the mode in the slits 
to the outgoing plane wave, that are computed using 
RCWA. While such an approach has proved extremely 
accurate in th e past for t he Extraordinary Optical Trans¬ 
mission (EOT)|iAJ is | 2i | 29| quite unexpectedly, it fails 
to predict the position of the resonance given by local 
RCWA full simulations (see Fi g. [^ . 

The assumption underlying ra is that only the prop¬ 
agative mode is able to reach the bottom of the grooves. 
However, because the depth of the grating is smaller than 
the skin depth, not all the evanescent modes are atten¬ 
uated enough to be neglected. To be more precise, one 
mode in particular, although it is not propagative in the 
z direction, presents an attenuation constant that is so 
low that it is still significantly strong at the bottom of 
the slits. 

The profile of this mode is shown Fig. [^and it is rela¬ 
tively flat. This mode is thus very close to the attenuated 
wave that is excited in the metal when a plane wave is re¬ 
flected by a metallic screen. Moreover, each reflection of 
the propagating mode at the top or at the bottom of the 
grooves generates this mode too, so that it fully partici¬ 
pates to the resonance. We propose here a very general 
two-modes approach to model the reflection coefficient of 
the structure. 

In this framework, we call Ai and Bi (resp. A 2 and 
B 2 ) are the amplitude of the propagating (resp. evanes¬ 
cent) mode traveling (resp. attenuating) downwards or 
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The downward amplitudes can be written as a result of 
the direct excitation by the incoming plane wave, in ad¬ 
dition to the reflection of the modes inside the slits 

r = toi + ri Bi + r2i B 2 

\ ^2 = to2 + ?'2 ^2 -I- ri2 -B 2 

and the equivalent is obtained at the bottom of the 
grooves for the upward amplitudes 

B^ = r\ Ai -f B 2 

B 2 24i -b rj2 ^2 

This forces us to introduce all the above new reflection 
and coupling coefficients, that can be computed using 
RCWA. This approach is inspired by mode recycling in 
photonic crystal cavitie^^, generalized here to (i) a non 
symmetrical plasmonic cavity and (ii) the case where the 
evanescent mode can be directly excited by the incident 
wave. 

We first use the relation above to eliminate B 2 from 
the whole system, yielding 



' r = ro + tioSi + t2o{rlA2e^^^^’^ + 

Ai = <01 + + r 2 i(r^ ^262*'=='' + 

A 2 = t 02 + ri(r^A2e2*'==^ + + ri 2 Bi 

. Bi = -b 

We finally get from these equations an expression for 
A 2 

A 2 = ^02 “t“ cAi + 

where 



Wavelength (nm) 

FIG. 2. (Color online) (a) Representation of the general¬ 
ized cavity model with the different reflection coefficients, (b) 
Top: Profile of the propagating mode (black curve) and the 
first evanescent mode(red curve). Bottom: Diverging effec¬ 
tive index of the propagating mode as a function of the gap 
width (solid line) and the mode guided in a single isolated slit, 
the gap-plasmon polariton (dotted line). (c) Reflection coeffi¬ 
cient of the whole structure, considering only the fundamental 
Bloch mode (black curve), a 2 modes model (red curve) and a 
complete simulation (green curve) for different groove depth 
(left: 7nm; center: 3 nm; right: 2 nm). Inset: simulation with 
Comsol multiphysics of the magnetic field and the Poynting 
vector at the entrance of the slits. 


,/ ^ ^02 

02 - 1 _ ^^j,b^2ik2h 

1 — 7'2r2e^*^2^ 

, _ ?^12 

12 2 ^ _ ^^^b^2ik2h 

We can thus replace A 2 in the above equations, to yield 

r = Tg -b t'^^Bi + aAi 
Ai = toi + f'lBi + i"2i(r2(tQ2 + cAi 

^ Bi = + r^,(t'o2 + cAi + 


upwards (see Fig. [^. The reflection coefficient of the 
whole structure can then be written as the result of the 
reflection on the metallic plane plus the light that comes 
from the two modes inside the grating layer 


where 


r = To -btio.Bi +^20 ^ 2 - 


(4) 


Tg = ro -b t20?’2^02e^*''"'' 
^10 = ^10 + ^20^2^12®^*^^^ 
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And finally, the system reduces to 

( r = rg + t'lgSi + aAi 
J Ai = r[Bi + Iqi 


where 


1 - 


toi = toi + t’o2f2irle^'’’"'"^ 

^b' _ ‘1^21^^ _ 


It is not obvious yet that these equations can be used 
to yield a generalized Fabry-Perot formula, because two 
terms have appeared that do not exist in the classical 
one-mode model. More precisely, to retrieve the exact 
same equations, tg 2 and a would have to vanish. This is 
usually the case in photonic crystals - but physically here, 
the evanescent mode is directly and efficiently excited by 
the incoming wave so that fg 2 is not negligible. 

We now replace Ai by its expression in the reflexion 
coefficient’s formula, to yield 

{ r = rg -f f'l'gBi 
Ai = r'lBi + tgi 


Then the reflection coefficient can be put under a form 
similar to the Fabry-Perot formula 


r = r, 


e// 


teffTi e 


b'2ik\h 


^ fj^b ^2'ik^h 


( 12 ) 


This formula can be considered as a generalized Fabry- 
Perot formula. The agreement between this formula and 
a full RCWA simulation (see Fig. (c)), that can be 
considered as a multi-mode model, is excellent. The two 
modes are thus the only ones that are responsible for 
the resonance. More precisely, the propagating mode is 
responsible for the resonance and the evanescent mode is 
responsible for a shift of this resonance compared to the 
one-mode model. 

The effective reflection coefficients are easy to com¬ 
pute, once the real coefficients are extracted from the 
scattering matrices of the interfaces between the differ¬ 
ent space region^^. The resonance condition now reads 


arg(rj) -|- arg(rj ) -I- 2^{ki)h = 2rmr (13) 

where m is a relative integer. Fig. shows the phase of 
the effective reflection coefficients, compared to the phase 
of the real reflection coefficients. There is essentially a 
shift of the phase of r[ compared to the one of ri. This 
totally explains why, compared to the one-mode model 
prediction, the resonance is shifted since the shift in the 
phase has a direct impact on the resonance condition. 
The fact that the slope of the phase is not changed means 
that there is almost no impact on the quality factor of the 
resonance, so that the only impact of the second mode 
on the resonance is to shift it without changing its width. 


where new effective coefficients are introduced: 


rg = r[) -b at'gi 

(5) 

^10 = ^10 + “^1 

( 6 ) 

Finally, Ai can be written 


1 _ ^01 + ^02^1 
^ ^ _ j,f^j.b' ^2ikih 

( 7 ) 

and the reflection coefficient, by eliminating Bi gives 

r = rg -b ti'grf -b t'/gC 

(8) 

leading to the desired result 


„ ,/ b' 2ikih^" 1^1 I /'I 

_" , 4." 4.” , C;^g^I^g -1- rg2r^J 

' 'o Go*^ 02 -r ■ ■ \ ^iiki'h 

( 9 ) 

If we call 


■'"e// = ■'’o + ^10^02 

( 10 ) 

^e// = ^lo(^10 + ^02^l) 

( 11 ) 
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FIG. 3. (Color online) Phase of the reflection coefficients. 
Top: Phase of ri and r[ (black and green curves respectively). 
Bottom: Phase of rj and r\ (black and green lines resp.). 


III. IMPACT OF SPATIAL DISPERSION 


Deeply subwavelength gratings are very interesting be¬ 
cause the absorbing resonance is the sign that a guided 













5 


mode with a very high wavevector has been efficiently 
excited in the structure. Such modes are likely to be 
influenced by the repulsion between electron^ because 
their effective wavelength is so small that is approaches 
the mean free path of electrons in the metaPI^, which is 
the relevant scale for nonlocality. That is the reason why 
the spatial dispersion is expected to have an impact on 

the slit array’s response. _ 

We use here the hydrodynamic model^ to take the 
intrinsic nonlocality of the metallic response into account. 
The currents J corresponding to the movement of the free 
electrons trapped in the metal can be taken into account 
as an effective polarization P/ defined ^ P/ = J. The 
electron gas can be considered as a fluicP, leading to the 
following linearized equation 

-/32V(V.P/)+ P/+7Py = eoW^E (14) 

where 7 is the damping factor, uip the plasma frequency 
and P ~ 1.35.10®m/s. All these parameters (except /3), as 
well as the dispersive susceptibility Xb due to interband 
transitions are taken from careful fits of the available 
experimental data using a Drude and Brendel-Bormann 
modeP^. This allows us to clearly distinguish between 
the response of the jellium and the response of the back¬ 
ground, that we assume is local. Maxwell’s equations 
then reduce tcP 


V X E = iw/roH 
V X H = -*a;eo(l + Xb)E + Pf 


where, thanks to (14) the polarization can be written as 
,2 


P/ = 


eo-w. 


P 




E-(1 + X6)4v(V-E) 


UJ. 


P 


(15) 


The resolution of such equations in a multilayered 
structure requires the introduction of an additional 
boundary condition (ABC). The most obvious ABC is 
in that case to consider that no electrons are allowed to 
get out of the metaP^, which means J.n = P/.n = 0, 
where n is the normal to the interface - and it turns out 
to be one of the most conservative ABC, so that nonlocal 
effects are not likely to be overestimatecP. 

In this framework the dispersion relation giving the 
propagation constant of the modes is modifiecP^ and be¬ 
comes 


1 - 
1 

-f - 


eQ sinh(K(e) 


kz sinh(K/e) 
Kt ena 
ena 1 — Kt 

sinh(K/e) 


= cosh(Ka) cosh(Kte) 
/ e 




\PKteJ _ 


sinh(Ka) sinh(Kje) 


sinh(Ka) 


P 

esinh(Kte) 


(1 — cosh(Kte) cosh(Kie) — 


Kt 


cosh(Ka) cosh(tt;te) 


(16) 


where e = d-a,n= §(i-Y^) , Kf = + 

^ , Kp = k1 — ek^ and k? = k^ — ko. This allows 


l+Xb 


to consider the impact of nonlocality on the propagation 
constants of the two modes that are involved in the res¬ 
onance of the structure - both the propagating and the 
least evanescent of the remaining modes, as shown Figure 
The figure actually shows that the propagating mode 
significantly more sensitive to nonlocality because of its 
large wavevector fci. A high wavevector actually means 
a large value for ft and thus a noticeable effect on the 
dispersion relation. 




FIG. 4. (Color Online) Impact of nonlocality on the effec¬ 
tive index (phase index) of the propagating mode (solid line) 
and on the evanescent mode (dashed line). The real (top) 
and imaginary (bottom) parts of the propagation constant 
are shown for a purely local approach (in black) and in the 
framework of the hydrodynamic model (red). 


In order to assess the impact of nonlocality on the 
whole structure, we have used COMSOL simulations, as 
a full nonlocal modal method is still beyond our reach. 
As expected, the resonance is blue-shifted compared to a 
local calculation (from 15 nm for 3 nm wide slits, see 
(b) to 24 nm for 2 nm wide slits, see Fig. (a)). Inter¬ 
estingly, using a two-mode model but changing only the 
propagation constant ki of the fundamental mode as com¬ 
puted using the nonlocal dispersion relation (161 allows 
to account for most of the shift (see Fig. [^. This demon¬ 
strates that the major reason why the whole structure is 
sensitive to the spatially dispersive response of the metal 
is that nonlocality has an impact on the wavevector of 
the single mode propagating in the slits. Its wavevector 
is actually not as high as predicted by the local theory, 
thus leading to a blue shift of the resonance. 


IV. CONCLUSION 

The mode propagating in a single extremely thin slit 
is often called a gap-plasmon polariton, to distinguish 
it from the surface plasmon polariton. The propagating 
mode considered here is of course closely related and is 
submitted to the same physical effect: the presence of the 
metal makes the mode slow, with a very large wavevector. 
This allows to understand (i) that deeply subwavelength 
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FIG. 5. (Color Online) Reflection coefficient computed using 
COMSOL (fully nonlocal calcnlation, blue curves), using a 
local 2-modes model (red curves) and using a nonlocal version 
of the 2-modes model (green curves) for 3 nm wide grooves 
(top) and 2 nm wide grooves (bottom). 


structures, sometimes smaller than the skin depth, can 
still be considered as cavities and (ii) why, as we have 
shown here, the smallest resonators are likely to be in¬ 
fluenced by spatial dispersion in metals. This class of 
resonators are called gap-plasmon resonators and it has 
recently been demonstrated experimentally that these 
resonators could constit ute extra ordinarily efficient con¬ 
centrators and absorber j^^b6 | 33H l5] a nd pr oduce totally 
unprecedented Purcell enhancementpSEII Theoretical 
studies show that they have a lo t of potential as scat- 
terers for light extraction tocPsMI. We are thus confident 
that the very general theoretical tools we have introduced 
here can be useful to study these structures and help de¬ 
sign them in the future, as they reach sizes that are well 
below the size of conventional cavity resonators. 
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